# Clear memory 
rm(list = ls())

# Load packages
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  data.table, janitor, magrittr, fixest,
  broom, tidyverse, tidylog, kableExtra, modelsummary
)
options("tidylog.display" = NULL)

# Import function to write table
source("scripts/4-empirical-results/1-0-create-latex-fragment-function.R")

#====================
## SET UP DATA
#====================

# map to rename modelsummary tables to match etable
gm <- tibble::tribble(
  ~raw,        ~clean,          ~fmt,
  "nobs",      "Observations",  0
)

amenity_df <- fread("data/amenity_panel.csv")
emissions_df <- fread("data/emissions.csv")

small_emissions_df = emissions_df %>% filter(pollutant %in% c("PM25", "VOC", "NH3", "SO2", "NOX"))

#====================
## EMISSIONS: MAIN
#====================

# Start with single
regs_single <- list()

regs_single[["(1)"]] <-
  fepois(
    amount ~ i(new_post_nonattainment, ref = 1) |
      fips + year + pollutant,
    data = small_emissions_df, 
    cluster = small_emissions_df$state
  )

regs_single[["(2)"]] <-
  fepois(
    amount ~ i(new_post_nonattainment, ref = 1)  |
      fips + year + pollutant^year,
    data = small_emissions_df, 
    cluster = small_emissions_df$state
  )

regs_single[["(3)"]] <-
  fepois(
    amount/annual_total_wages ~ i(new_post_nonattainment, ref = 1)  |
      fips + year + pollutant,
    data = small_emissions_df, 
    cluster = small_emissions_df$state
  )

regs_single[["(4)"]] <-
  fepois(
    amount/annual_total_wages ~ i(new_post_nonattainment, ref = 1) |
      fips + pollutant^year,
    data = small_emissions_df, 
    cluster = small_emissions_df$state
  )

# Create model dictionary that can be used for labeling variables.
dictionary <- c(
  "new_post_nonattainment = 0" = "\\addlinespace \\hspace{.5cm} $\\beta^p_{\\eta}$",
  new_post_nonattainment = "\\addlinespace \\hspace{.5cm} $\\beta^p_{\\eta}$",
  fips = "County",
  year = "Year",
  pollutant = "Pollutant"
)

# Top
create_latex_fragment(
  etable_model_list = regs_single,
  variable_dictionary = dictionary,
  file_name = "output/table_1_emissions_table_panel.tex",
  group_header = "\\emph{A. Combined}",
  position = "top",
  keep_list = "beta",
  append_existing = FALSE,
  remove_text = "\\$=\\$ 0"
)

# Now do each pollutant

regs <- list()

regs[["(1)"]] <-
  fepois(
    amount ~ i(new_post_nonattainment, pollutant, ref = 1) |
      fips + year + pollutant,
    data = small_emissions_df, 
    cluster = small_emissions_df$state
  )

regs[["(2)"]] <-
  fepois(
    amount ~ i(new_post_nonattainment, pollutant, ref = 1)  |
      fips + year + pollutant^year,
    data = small_emissions_df, 
    cluster = small_emissions_df$state
  )

regs[["(3)"]] <-
  fepois(
    amount/annual_total_wages ~ i(new_post_nonattainment, pollutant, ref = 1) |
      fips + year + pollutant,
    data = small_emissions_df, 
    cluster = small_emissions_df$state
  )

regs[["(4)"]] <-
  fepois(
    amount/annual_total_wages ~ i(new_post_nonattainment, pollutant, ref = 1) |
      fips + pollutant^year,
    data = small_emissions_df, 
    cluster = small_emissions_df$state
  )

# Export coefs for model
coefs = regs[["(4)"]]$coefficients
coefs = data.frame(t(coefs))
fwrite(coefs, file = "data/counterfactual/inputs/emissions_reductions.csv") 

# Create model dictionary that can be used for labeling variables.
dictionary <- c(
  "new_post_nonattainment = 0" = " ",
  NH3 = "Ammonia ($\\beta^{NH_3}_{\\eta}$)",
  NOX = "Nitrogen Oxides ($\\beta^{NO_x}_{\\eta}$)",
  PM25 = "Fine Particulates ($\\beta^{PM_{2.5}}_{\\eta}$)",
  SO2 = "Sulfur Dioxide ($\\beta^{SO_2}_{\\eta}$)",
  VOC = "Volatile Organics ($\\beta^{VOC}_{\\eta}$)",
  new_post_nonattainment = "\\addlinespace \\hspace{.5cm}",
  fips = "County",
  year = "Year",
  pollutant = "Pollutant"
)

# Show results in base etable
etable(regs, dict = dictionary, keep = c("beta", "hspace"), interaction.combine = "~")

# Bottom
create_latex_fragment(
  etable_model_list = regs,
  variable_dictionary = dictionary,
  file_name = "output/table_1_emissions_table_panel.tex",
  group_header = "\\emph{B. By Emitted Pollutant}",
  position = "bottom",
  keep_list = "beta",
  append_existing = TRUE,
  remove_text = "\\$=\\$ 0~Pollutant \\$=\\$",
  extra_text = list("-Outcome" = c("\\multicolumn{1}{c}{Emissions}", "\\multicolumn{1}{c}{Emissions}", "\\multicolumn{1}{c}{Emissions Intensity}", "\\multicolumn{1}{c}{Emissions Intensity}"))
)

#====================
## EMISSIONS: ROBUST
#====================

regs_single <- list()

small_emissions_df_rob = small_emissions_df %>% filter(year <= 1996)

regs_single[["(1)"]] <-
  fepois(
    amount/annual_total_wages ~ i(new_post_nonattainment, ref = 1) |
      fips + pollutant^year,
    data = small_emissions_df_rob,
    cluster = small_emissions_df_rob$state
  )

small_emissions_df_rob = small_emissions_df %>% filter(year <= 1997)

regs_single[["(2)"]] <-
  fepois(
    amount/annual_total_wages ~ i(new_post_nonattainment, ref = 1) |
      fips + pollutant^year,
    data = small_emissions_df_rob,
    cluster = small_emissions_df_rob$state
  )

small_emissions_df_rob = small_emissions_df %>% filter(year <= 1998)

regs_single[["(3)"]] <-
  fepois(
    amount ~ i(new_post_nonattainment, ref = 1) |
      fips + pollutant^year,
    data = small_emissions_df_rob,
    cluster = small_emissions_df_rob$state
  )

small_emissions_df_rob = small_emissions_df %>% filter(year <= 1999)

regs_single[["(4)"]] <-
  fepois(
    amount/annual_total_wages ~ i(new_post_nonattainment, ref = 1) |
      fips + pollutant^year,
    data = small_emissions_df_rob,
    cluster = small_emissions_df_rob$state
  )

small_emissions_df_rob = small_emissions_df %>% filter(year <= 2000)

regs_single[["(5)"]] <-
  fepois(
    amount/annual_total_wages ~ i(new_post_nonattainment, ref = 1) |
      fips + year + pollutant,
    data = small_emissions_df_rob,
    cluster = small_emissions_df_rob$state
  )

small_emissions_df_rob = small_emissions_df %>% filter(year <= 2001)

regs_single[["(6)"]] <-
  fepois(
    amount/annual_total_wages ~ i(new_post_nonattainment, ref = 1) |
      fips + pollutant^year,
    data = small_emissions_df_rob,
    cluster = small_emissions_df_rob$state
  )

rows <- tribble(
  ~term, ~"(1)", ~"(2)", ~"(3)", ~"(4)", ~"(5)", ~"(6)",
  "County FE", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes",
  "Year FE", "No", "No", "No", "No", "No", "No",
  "Pollutant FE", "No", "No", "No", "No", "No", "No",
  "Pollutant-Year FE", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes",
  "Latest Year", "1996", "1997", "1998", "1999", "2000", "2001"
)


msummary(regs_single,
  cluster = c("fips"),
  coef_map = c(
    "new_post_nonattainment::0" = "$\\beta^p_{\\eta}$"
  ),
  stars = c("*" = .1, "**" = .05, "***" = .01),
  gof_omit = "R2|AIC|BIC|Log.|FE|Std",
  gof_map = gm,
  fmt = "%.2f",
  add_rows = rows,
  escape = FALSE,
  title = "Difference-in-differences estimates of the effect of nonattainment on the implicit emissions price varying the sample period.\\label{tab:emissions_length}",
  notes = list("The coefficients are estimated using data from 1990, and from 1996 to the year listed in the Latest Year row. The 1991--1995 gap reflects years in which NEI data are not available. Robust standard errors are clustered at the state level."),
  output = "output/table_d1_emissions_table_length.tex"
)

#====================
## AMENITIES
#====================

regs <- list()

regs[["(1)"]] <-
  fepois((share / initial_own_share) ~
  I(as.numeric(treated_dest) - as.numeric(treated_orig)) |
    origin^destination + year, amenity_df,
  cluster = ~ origin_st + destination_st
  )

regs[["(2)"]] <-
  fepois((share / initial_own_share)/((dest_wage / origin_wage) / (dest_ma / origin_ma)^(-1 / 4)) ~
  I(as.numeric(treated_dest) - as.numeric(treated_orig)) |
    origin^destination + year, amenity_df,
  cluster = ~ origin_st + destination_st
  )

regs[["(3)"]] <-
  fepois((share / initial_own_share)/((dest_wage / origin_wage) / (dest_ma / origin_ma)^(-1 / 4)) ~
  treated_dest + treated_orig |
    origin^destination + year, amenity_df,
  cluster = ~ origin_st + destination_st
  )

rows <- tribble(
  ~term, ~"(1)", ~"(2)", ~"(3)", 
  "Origin-Destination FE", "Yes", "Yes", "Yes", 
  "Year FE", "Yes", "Yes", "Yes", 
  "Real Wage", "Omitted", "Coef. Fixed", "Coef. Fixed",
  "Origin vs Destination", "Fixed", "Fixed", "Free"
)


msummary(regs,
  cluster = c("origin", "destination"),
  coef_map = c(
    "I(as.numeric(treated_dest) - as.numeric(treated_orig))" = "1(Dest. Nonattainment) - 1(Orig. Nonattainment)",
    "treated_destTRUE" = "1(Destination Nonattainment)",
    "treated_origTRUE" = "1(Origin Nonattainment)"
  ),
  stars = c("*" = .1, "**" = .05, "***" = .01),
  gof_omit = "R2|AIC|BIC|Log.|FE|Std",
  gof_map = gm,
  fmt = "%.3f",
  add_rows = rows,
  title = "Difference-in-differences estimates of the effect of nonattainment on local amenities.\\label{tab:amenities}",
  notes = list("Robust standard errors are clustered two ways at the origin state and destination state levels."),
  output = "output/table_b1_amenity_table.tex"
)